t0 <- Sys.time()
set.seed(2) # make workflow repeatable
library(tidyverse)
library(ggplot2)
library(ggExtra)
library(ggpubr)
library(kableExtra)
library(tidyr)
library(readxl)
library(lubridate)
library(sf)
library(raster)
library(amt)
library(purrr)
library(runjags)
library(patchwork)
library(lme4)
library(lmerTest)
library(emmeans)
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(DHARMa)
select<-dplyr::select
crs_utm18n <- "EPSG:26918" # UTM zone 18N, https://epsg.io/32618
lu <- function(x) length(unique(x)) # a helper function to save space
# fish metadata recorded after tagging and release
fishids <- read_excel("BowfinInventory_gj.xls") %>% 
  mutate(id=as.character(FishID))
# fish observation data
fishobs_tab <- read_csv("bowfin_data.csv") %>% 
  left_join(fishids) %>% 
  filter(!is.na(Latitude)) %>% # filter out missing coordinates
  mutate(DOY=yday(DateTime)) %>% # need DOY of all datetime observations
  filter(!is.na(distday)&Distance>0) %>% # keep only observations of movement
  mutate(Sex=factor(Sex), Year=factor(Year), ReleaseSite=factor(ReleaseSite),
         Zdoy=scale(.$endDOY)[,1])
fishobs <- fishobs_tab %>% filter(!is.na(Latitude)) %>% 
  st_as_sf(coords=c("Longitude", "Latitude"), crs=4326) %>%
  st_transform(crs=crs_utm18n)
# Bayesian change-point models and results
load(file="mcp_out.RData")
# Bayesian change-point results summery table
mcp_summary <- read_excel("MCPresultstab_20240303.xlsx") %>% 
  select(-Notes) %>% 
  left_join(fishids)
# sf object of initial capture locations for each fish
caploc <- fishids %>%
  mutate(x=if_else(ReleaseSite=="Bradburys", -76.11321, -75.92999),
         y=if_else(ReleaseSite=="Bradburys", 43.25366, 43.17359)) %>%
  st_as_sf(coords=c("x", "y"), crs=4326) %>%
  st_transform(crs=crs_utm18n) %>%
  mutate(Year=year(ReleaseDate), DOY=yday(ReleaseDate)) %>%
  select(FishID, Year, DOY) 

Supporting spatial data

Load spatial datasets and project everything into NAD83/UTM zone 18N.

# Oneida Lake boundary from NHD v2.1
olb <- sf::st_read(dsn="supporting", layer="OneidaLake_utm18n") 
## Reading layer `OneidaLake_utm18n' from data source 
##   `C:\Users\gj93\Box\Bowfin_telemetry\bowfin_analysis_mobile\bowfin_HR_SignerWorkflow_220709\manuscript_version\supporting' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 14 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 407400 ymin: 4777457 xmax: 440730.1 ymax: 4789914
## Projected CRS: NAD83 / UTM zone 18N
lrwpoly <- sf::st_read(dsn="supporting",
                       layer="lrw_poly",
                       crs=crs_utm18n) 
## Reading layer `lrw_poly' from data source 
##   `C:\Users\gj93\Box\Bowfin_telemetry\bowfin_analysis_mobile\bowfin_HR_SignerWorkflow_220709\manuscript_version\supporting' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 400003.9 ymin: 4773813 xmax: 444803.9 ymax: 4793563
## Projected CRS: NAD83 / UTM zone 18N
bathy <- raster("supporting/depth_2")
rw <- raster("supporting/rwindicator")
lk <- raster("supporting/oneidaind")
lrw <- raster("supporting/lrwindicator")
# create data frames from rasters, clipped to lake shore, for plotting
bathy_df <- as.data.frame(bathy*lk, xy=TRUE) 
names(bathy_df)[3] <- "value"
rw_df <- as.data.frame(rw*is.na(lk), xy=TRUE) %>% 
  mutate(value=ifelse(layer==1, 1, NA))
lrw_df <- as.data.frame(lrw*!is.na(lrw), xy=TRUE) %>% 
  mutate(value=ifelse(layer==1, 1, NA))
lrw_mask_df <- mutate(lrw_df, value=ifelse(is.na(value),1,NA))

Map locations across available habitats in Oneida Lake and surrounding wetland and stream habitats.

# plot locations by FishID - base plot for fig 1 and results figures
plotlocations <- ggplot() + 
  geom_tile(data=rw_df %>% filter(!is.na(value)), 
            aes(x=x, y=y), fill='gray90', color='gray90') +
  geom_tile(data=bathy_df %>% filter(!is.na(value)), 
            aes(x=x, y=y, fill=value, color=value)) +
  scale_fill_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
  scale_color_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
    geom_sf(data=fishobs, shape=3, alpha=0.5) +
  geom_sf(data=caploc, shape=19, col="gold") +
  coord_sf(xlim=c(408000, 440000), ylim=c(4773813, 4793563)) +
    theme_void() + #theme(legend.position="none")+
    ggspatial::annotation_north_arrow(which_north="true", pad_y=unit(0.25, "in"),
        style = ggspatial::north_arrow_fancy_orienteering) +
    ggspatial::annotation_scale(location = "bl", width_hint = 0.5) + 
  theme(plot.title=element_text(hjust=0.5))

Depth and Distance from shore.

Visualize and analyze sex differences.

fishobs_tab %>% 
  group_by(FishID, Sex) %>% 
  summarize(ShoreDist_m=mean(ShoreDist_m, na.rm=TRUE),
            Depth_m=mean(Depth_m, na.rm=TRUE)) %>% 
  ungroup() %>% 
  group_by(Sex) %>% 
  summarize(mean_ShoreDist_m=mean(ShoreDist_m, na.rm=TRUE),
            sd_ShoreDist_m=sd(ShoreDist_m, na.rm=TRUE),
            n_ShoreDist_m=sum(!is.na(ShoreDist_m)),
            mean_Depth_m=mean(Depth_m, na.rm=TRUE),
            sd_Depth_m=sd(Depth_m, na.rm=TRUE),
            n_Depth_m=sum(!is.na(Depth_m)))
## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.
## # A tibble: 2 × 7
##   Sex   mean_ShoreDist_m sd_ShoreDist_m n_ShoreDist_m mean_Depth_m sd_Depth_m
##   <fct>            <dbl>          <dbl>         <int>        <dbl>      <dbl>
## 1 F                 130.           48.1            18         1.58      0.371
## 2 M                 109.           63.0            18         1.35      0.227
## # ℹ 1 more variable: n_Depth_m <int>

Repeated measures analysis of depth by sex.

lmer(log(Depth_m) ~ Sex + (1|FishID), data=fishobs_tab) %>% 
  summary(ddf="Satterthwaite")
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(Depth_m) ~ Sex + (1 | FishID)
##    Data: fishobs_tab
## 
## REML criterion at convergence: 1466.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6265 -0.6147 -0.0067  0.6471  2.6218 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FishID   (Intercept) 0.01494  0.1222  
##  Residual             0.30774  0.5547  
## Number of obs: 864, groups:  FishID, 36
## 
## Fixed effects:
##             Estimate Std. Error       df t value Pr(>|t|)    
## (Intercept)  0.28729    0.04372 34.68911   6.570 1.44e-07 ***
## SexM        -0.14765    0.05857 31.17784  -2.521    0.017 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## SexM -0.747
lmer(log(Depth_m) ~ Sex + (1|FishID), data=fishobs_tab) %>% 
  anova(ddf="Satterthwaite")
## Type III Analysis of Variance Table with Satterthwaite's method
##     Sum Sq Mean Sq NumDF  DenDF F value  Pr(>F)  
## Sex 1.9558  1.9558     1 31.178  6.3553 0.01703 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Repeated measures analysis of distance from shore by sex.

lmer(log(ShoreDist_m+1)  ~ Sex + (1|FishID), data=fishobs_tab) %>% 
  summary(ddf="Satterthwaite")
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log(ShoreDist_m + 1) ~ Sex + (1 | FishID)
##    Data: fishobs_tab
## 
## REML criterion at convergence: 2686.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.7910 -0.4359  0.1155  0.6736  2.2524 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  FishID   (Intercept) 0.1727   0.4155  
##  Residual             1.4405   1.2002  
## Number of obs: 824, groups:  FishID, 36
## 
## Fixed effects:
##             Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)   4.4551     0.1251 33.9960  35.619  < 2e-16 ***
## SexM         -0.4843     0.1698 31.3363  -2.853  0.00761 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##      (Intr)
## SexM -0.737
lmer(log(ShoreDist_m+1)  ~ Sex + (1|FishID), data=fishobs_tab) %>% 
  anova(ddf="Satterthwaite")
## Type III Analysis of Variance Table with Satterthwaite's method
##     Sum Sq Mean Sq NumDF  DenDF F value   Pr(>F)   
## Sex 11.722  11.722     1 31.336  8.1375 0.007614 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
depthdistmvmntplot <- fishobs_tab %>% 
  group_by(FishID, Sex) %>% 
  summarize(barShoreDist=mean(ShoreDist_m+1, na.rm=TRUE),
            barDepth_m=mean(Depth_m, na.rm=TRUE),
            barDistDay=mean(distday, na.rm=TRUE))
## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.
pShoreDist <- depthdistmvmntplot %>% 
  ggplot(aes(barShoreDist, color=factor(Sex))) +
  scale_color_grey() +
  geom_density(size=1) +  
  scale_x_log10() +
  ylab("Density") +
  xlab("Shore Distance (m)") + 
  theme_classic() + theme(legend.position = "none")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
pDepth <- depthdistmvmntplot %>% 
  ggplot(aes(barDepth_m, color=factor(Sex))) +
  scale_color_grey() +
  geom_density(size=1) + 
  scale_x_log10() +
  ylab("") +
  xlab("Depth (m)") + 
  theme_classic() + theme(legend.position = "none")
pdistday <- depthdistmvmntplot %>% 
  ggplot(aes(barDistDay, color=factor(Sex))) +
  scale_color_grey() +
  geom_density(size=1) + 
  scale_x_log10(breaks=c(1,10,100,1000)) +
  ylab("") +
  xlab("Movement (m/day)") + 
  theme_classic() + labs(color="Sex") 
fig3 <- pShoreDist + pDepth + pdistday + 
  plot_annotation(tag_levels=list(c("a)", "b)", "c)"))) +
  plot_layout(widths=c(1,1,1), ncol=3, nrow=1)  
ggsave("fig3_freq.png", fig3, width=6.5, height=2, units="in", dpi=600, 
       bg="white")
knitr::include_graphics("fig3_freq.png")

Movement analysis

Analysis of minimum distance moved per day by sex, year, and day of year.

# population average and sd of individual average movement rate
mMdat <- filter(fishobs_tab, dayssince<22)
group_by(mMdat, FishID, Sex) %>% 
  summarize(meandistday=mean(distday, na.rm=TRUE)) %>% 
  group_by(Sex) %>% 
  summarize(meanM=mean(meandistday),
            sdM=sd(meandistday))
## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.
## # A tibble: 2 × 3
##   Sex   meanM   sdM
##   <fct> <dbl> <dbl>
## 1 F      221.  157.
## 2 M      114.  107.
# population average of within-individual SD of movement rate
mMdat <- filter(fishobs_tab, dayssince<22)
group_by(mMdat, FishID, Sex) %>% 
  summarize(sddistday=sd(distday, na.rm=TRUE)) %>% 
  group_by(Sex) %>% 
  summarize(meansdM=mean(sddistday, na.rm=TRUE))
## `summarise()` has grouped output by 'FishID'. You can override using the
## `.groups` argument.
## # A tibble: 2 × 2
##   Sex   meansdM
##   <fct>   <dbl>
## 1 F        280.
## 2 M        253.
mM <- lmer(log(distday) ~ Sex*Year + Sex*Zdoy + (Zdoy|FishID), data=mMdat)
anova(mM, ddf="Satterthwaite")
## Type III Analysis of Variance Table with Satterthwaite's method
##          Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)    
## Sex      64.961  64.961     1  34.91 21.3413 5.057e-05 ***
## Year      7.585   3.793     2 744.60  1.2459  0.288272    
## Zdoy      8.989   8.989     1  31.25  2.9531  0.095606 .  
## Sex:Year 30.519  15.260     2 744.60  5.0131  0.006876 ** 
## Sex:Zdoy 12.960  12.960     1  31.25  4.2576  0.047465 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
margMdf <- emmip(mM, ~ Sex|Year|Zdoy, lmer.df="satterthwaite", plotit=FALSE)
pd <- position_dodge(.3)
p <- ggplot(data=margMdf, aes(x=Year, y=exp(yvar), fill=Sex, Sroup=Sex)) +
  geom_errorbar(aes(ymin=exp(yvar-SE), ymax=exp(yvar+SE)), 
                color='black', width=.2, position=pd) +
  geom_line(position=pd) +
  geom_point(color='black', size=8, position=pd, shape=21) +
  scale_fill_grey() +
  ylab(expression("Movement "~(m/day))) +
  theme_classic() 
ggsave("fig4_Mmeans.png", p, width=6.5, height=4, units="in", dpi=600, 
       bg="white")
## `geom_line()`: Each group consists of only one observation.
## ℹ Do you need to adjust the group aesthetic?
knitr::include_graphics("fig4_Mmeans.png")

tab_model(mM, df.method = "satterthwaite")
  log(distday)
Predictors Estimates CI p
(Intercept) 4.60 3.92 – 5.28 <0.001
Sex [M] -2.09 -2.98 – -1.20 <0.001
Year [2010] -0.60 -1.27 – 0.07 0.078
Year [2011] -1.02 -1.73 – -0.31 0.005
Zdoy -0.38 -0.70 – -0.07 0.017
Sex [M] × Year [2010] 0.99 0.13 – 1.86 0.024
Sex [M] × Year [2011] 1.44 0.54 – 2.35 0.002
Sex [M] × Zdoy 0.42 0.00 – 0.83 0.047
Random Effects
σ2 3.04
τ00 FishID 0.38
τ11 FishID.Zdoy 0.16
ρ01 FishID 0.44
ICC 0.15
N FishID 36
Observations 813
Marginal R2 / Conditional R2 0.095 / 0.231
simulationOutput <- simulateResiduals(fittedModel=mM, plot=FALSE)
testTemporalAutocorrelation(simulationOutput, time=mMdat$datetime)
## Warning: Unknown or uninitialised column: `datetime`.
## DHARMa::testTemporalAutocorrelation - no time argument provided, using random times for each data point

## 
##  Durbin-Watson test
## 
## data:  simulationOutput$scaledResiduals ~ 1
## DW = 1.9293, p-value = 0.3131
## alternative hypothesis: true autocorrelation is not 0
mMreduced <- lmer(log(distday) ~ Sex*Year + Sex*Zdoy + (1|FishID), data=mMdat)
anova(mMreduced, mM)
## refitting model(s) with ML (instead of REML)
## Data: mMdat
## Models:
## mMreduced: log(distday) ~ Sex * Year + Sex * Zdoy + (1 | FishID)
## mM: log(distday) ~ Sex * Year + Sex * Zdoy + (Zdoy | FishID)
##           npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)   
## mMreduced   10 3301.2 3348.2 -1640.6   3281.2                        
## mM          12 3292.8 3349.2 -1634.4   3268.8 12.416  2   0.002013 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Change point models

Change point model outputs are produced in a separate script, and summaries of their results are uploaded here for presentation and for use in the home range analysis.

mcp_summary %>% arrange(Usedmod)
## # A tibble: 36 × 18
##        i FishID Topmod Usedmod     n class      Frequency Sex    TLmm   Wkg
##    <dbl>  <dbl>  <dbl>   <dbl> <dbl> <chr>          <dbl> <chr> <dbl> <dbl>
##  1     6     20      1       0     5 Stationary      150. F       704  3.1 
##  2     9     23      2       0    15 Stationary      150. M       576  1.5 
##  3    28     43      3       0    34 Stationary      150. M       616  1.95
##  4    29     44      0       0    31 Stationary      150. M       581  1.75
##  5    31     47      0       0    33 Stationary      150. M       570  1.9 
##  6    34     32      0       0     8 Stationary      150. F       705  3.35
##  7    13     28      1       1    33 Migration       150. F       676  2.4 
##  8     1     12      2       2    24 Migration       150. M       545  1.38
##  9     5     19      2       2     7 Migration       150. F       674  3   
## 10     7     21      2       2    27 Migration       150. F       646  1.95
## # ℹ 26 more rows
## # ℹ 8 more variables: ReleaseDate <dttm>, ReleaseSite <chr>,
## #   SurgeryTime.min <chr>, Surgeon <chr>, BatteryEnd <dttm>, Omit <dbl>,
## #   Notes <chr>, id <chr>

Reorganize fixed effect coefficients from change point models.

# concatenate seasonal change points and filter out non-seasonal models
mcp_fixefs <- lapply(1:length(mcp_topmod), function(x){
  summary(mcp_topmod[[x]]$mcmc_post)$quantiles %>% data.frame() %>% 
    rename(median=X50.) %>% 
    mutate(name=row.names(.),
           id=names(mcp_topmod[x])) %>% 
    filter(!grepl("\\[", name))
}) %>% bind_rows() 
# extract random effects - median change point by year by fish
mcp_cp1yr <- lapply(1:length(mcp_topmod), function(x){
  summary(mcp_topmod[[x]]$mcmc_post)$quantiles %>% data.frame() %>% 
    rename(median=X50.) %>% 
    mutate(name=row.names(.),
           id=names(mcp_topmod[x])) %>% 
    separate("name", into=c("name", "year"), sep="\\[|\\]", extra="drop") %>% 
    filter(name=="cp_1"|name=="cp_1_year") %>% 
    select(id, name, year, median) %>% 
    mutate(cp1day=median+pull(filter(., name=="cp_1"), median)) %>% 
    filter(!is.na(year))
}) %>% bind_rows() 
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 5, 6, 7,
## 8].
## Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 5, 6, 7,
## 8].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 5, 6, 7,
## 8].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1, 2].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1, 2].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 2 rows [1, 2].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 6 rows [1, 2, 6, 7, 8,
## 9].
## Warning: Expected 2 pieces. Missing pieces filled with `NA` in 9 rows [1, 2, 6, 7, 11,
## 12, 13, 14, 15].
mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_1"& class=="Migration") %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median), 
            n=n())
## Joining with `by = join_by(id)`
##       mean       sd      min      max  n
## 1 162.1114 20.92017 118.5156 193.3826 30

How many fish with change point models supported, and do their change point dates differ by sex? The top models for some fish contained change point estimates but were rejected as still not supporting migratory behavior due to lack of fit. Because of this, we need cp_1 parameters from fish tracks classified as “Migration”.

mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_1"& class=="Migration") %>% 
  group_by(Sex) %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median), 
            n=n())
## Joining with `by = join_by(id)`
## # A tibble: 2 × 6
##   Sex    mean    sd   min   max     n
##   <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F      160.  20.8  128.  193.    16
## 2 M      164.  21.6  119.  193.    14
mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_1"& class=="Migration") %>% 
  {t.test(.$median ~ .$Sex, var.equal = FALSE)}
## Joining with `by = join_by(id)`
## 
##  Welch Two Sample t-test
## 
## data:  .$median by .$Sex
## t = -0.52384, df = 27.19, p-value = 0.6046
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -20.01236  11.86997
## sample estimates:
## mean in group F mean in group M 
##        160.2116        164.2827

How many of these fish had two change points between (to and from) two intervals? We want to compare cp_2; just extracting that parameter is enough to filter to the target subset.

mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_2") %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
            n=n())
## Joining with `by = join_by(id)`
##       mean       sd      min      max  n
## 1 263.5515 20.07628 211.2955 289.8022 19
mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_2") %>% 
  group_by(Sex) %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median), 
            n=n())
## Joining with `by = join_by(id)`
## # A tibble: 2 × 6
##   Sex    mean    sd   min   max     n
##   <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F      254.  22.4  211.  280.     9
## 2 M      272.  13.8  244.  290.    10
mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_2") %>% 
  {t.test(.$median ~ .$Sex, var.equal = FALSE)}
## Joining with `by = join_by(id)`
## 
##  Welch Two Sample t-test
## 
## data:  .$median by .$Sex
## t = -2.0673, df = 13.032, p-value = 0.05917
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -36.5923290   0.8009524
## sample estimates:
## mean in group F mean in group M 
##        254.1327        272.0284

How many of the “seasonal” fish tracks had variance-only changes points. We can look at this by extracting and summarizing the number of fish for which Usedmod equals 1.

mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_1"&Usedmod==1) %>% 
  group_by(Sex) %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
            n=n())
## Joining with `by = join_by(id)`
## # A tibble: 1 × 6
##   Sex    mean    sd   min   max     n
##   <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F      183.    NA  183.  183.     1

Here are results for fish with single change point models.

mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_1"&Usedmod==2) %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
            n=n())
## Joining with `by = join_by(id)`
##       mean       sd      min      max  n
## 1 158.9763 24.76114 118.5156 193.0988 10
mcp_fixefs %>% left_join(mcp_summary) %>% 
  filter(name=="cp_1"&Usedmod==2) %>% 
  group_by(Sex) %>% 
  summarize(mean=mean(median), sd=sd(median), min=min(median), max=max(median),
            n=n())
## Joining with `by = join_by(id)`
## # A tibble: 2 × 6
##   Sex    mean    sd   min   max     n
##   <chr> <dbl> <dbl> <dbl> <dbl> <int>
## 1 F      163.  23.2  133.  193.     5
## 2 M      155.  28.3  119.  188.     5

Seasons

We delineated season using the change point analysis and used these delineations to classify fish locations.

mcp_sf_lines <- lapply(mcp_sf2, function(x) {
  summarize(x, do_union=FALSE) %>% st_cast("LINESTRING")
  })
pview <- c(407400.0, 4777456.9, 427398.1, 4789913.8 )
p2fun <- function(i){
  name <- paste(names(mcp_sf2)[[i]])
  sex <- filter(fishids, FishID==as.numeric(names(mcp_sf2)[[i]])) %>% pull(Sex)
  tl <- filter(fishids, FishID==as.numeric(names(mcp_sf2)[[i]])) %>% pull(TLmm)
  p1 <- mcp_sf2[[i]] %>%
    ggplot() + 
    geom_tile(data=filter(lrw_df, !is.na(value)), aes(x=x, y=y), 
              fill='gray90') +
      geom_sf(data=olb, color="light blue", fill=NA) +
    geom_sf(data=mcp_sf_lines[[i]], color='gray70') +
    geom_sf(data=slice(mcp_sf2[[i]], which(t_==max(t_)|t_==min(t_))), 
            shape=c(1, 2), size=5) + 
    geom_sf(aes(color=doy, group=year)) + 
    scale_color_viridis_c(name="DOY") +
    theme_void() + 
    coord_sf(xlim=pview[c(1,3)], ylim=pview[c(2,4)]) 
  p2 <- mcp_sf2[[i]]  %>% 
    mutate(season = if (exists('season', where=.)) season else 5) %>% 
    ggplot(aes(x=t_, y=nsd)) +
    ggtitle("") + 
    geom_vline(xintercept=as.POSIXct(c("2009-01-01 00:00:00", 
                                       "2010-01-01 00:00:00",
                                       "2011-01-01 00:00:00",
                                       "2012-01-01 00:00:00")),
               lty=2, color="dark gray") +
    geom_line(aes(group=factor(year(t_)))) +
    geom_point(aes(fill=factor(season)), color="grey50", shape=21, size=3) + 
    scale_fill_grey(name="Season", na.value="white") +
    #scale_color_discrete(name="Year") + 
    ylim(0, 350) +
    coord_cartesian(xlim=as.POSIXct(c("2009-01-01 00:00:00", 
                                      "2012-01-01 00:00:00"))) +
    ylab(expression(NSD~(km^2))) + xlab("Time") +
    theme_classic() +
    theme(legend.position="none")
  p2marg <- ggMarginal(p2, type="histogram", margin="y")
  p <- (p1+p2marg)
  return(p)
}
p.seasonalM <- p2fun(12) #FishID #27 fit2, seasonal male
p.seasonalF <- p2fun(23) #FishID #38 fit3, seasonal female
p.seasonalVar <- p2fun(13) #FishID #28 fit1, var-seasonal female
p.aseasonal <- p2fun(29) #FishID #44 Stationary male
fig5 <- p.seasonalM[[1]] + p.seasonalM[[2]] +
  p.seasonalF[[1]] + p.seasonalF[[2]] + 
  p.seasonalVar[[1]] +  p.seasonalVar[[2]] +
  p.aseasonal[[1]] +  p.aseasonal[[2]] + 
  plot_layout(nrow=4, ncol=2) +
  plot_annotation(tag_levels='a', tag_suffix = ')') 
ggsave("fig5_pCPM.png", fig5, width=8, height=8, units="in", dpi=600, bg="white")
knitr::include_graphics("fig5_pCPM.png")

pcpdist <- filter(mcp_fixefs, name=="cp_1") %>% left_join(fishids) %>% 
  ggplot(aes(factor(Sex), median, color=TLmm)) + 
  geom_violin(fill="white") + 
  ylab("Change point (day of year)") +
  xlab("Sex") + 
  ggbeeswarm::geom_beeswarm(size=3) + 
  scale_color_viridis_c() + theme_minimal() + labs(fill=c("P(DiffInt)"))
## Joining with `by = join_by(id)`
pcpdist
## Warning: The following aesthetics were dropped during statistical transformation: colour
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
##   the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
##   variable into a factor?

Define separate tracks for separate seasons.

# all relocations
dat_hr <- mcp_sf2_tab %>% left_join(fishids) %>% left_join(mcp_summary) %>% 
  mutate(x_=st_coordinates(.)[,1], y_=st_coordinates(.)[,2]) %>% 
  st_drop_geometry() %>% as_tibble() %>%
  select(x_, y_, t_, sex=Sex, tl=TLmm, kg=Wkg, id, rel=ReleaseSite, 
         t_rel=ReleaseDate, surgeon=Surgeon, year, season) 
## Joining with `by = join_by(id)`
## Joining with `by = join_by(id, FishID, Frequency, Sex, TLmm, Wkg, ReleaseDate,
## ReleaseSite, SurgeryTime.min, Surgeon, BatteryEnd, Omit, Notes)`
# all seasons
dat_all <- select(dat_hr, -season) %>% 
  make_track(x_, y_, t_, id, sex, tl, kg, rel, crs=crs_utm18n) %>% 
  nest(data =c(x_, y_, t_))
# spring & fall seasons separately
dat_season <- filter(dat_hr, season%in%c(0,1)) %>% 
  make_track(x_, y_, t_, id, sex, tl, kg, rel, season, crs=crs_utm18n) %>% 
  nest(data =c(x_, y_, t_))

Home Ranges

Organize relocation data into tidy framework employed by package amt. To do so, condition into a movement track list - a tibble with one row per fish and a list data column containing relocations. Next, add a new list column that contains the home-range estimate for each animal. Estimate home ranges using 4 methods as Signer et al. (2015) did, truncating the dataset to only tracks with greater than 9 relocations.

hrlevel <- 0.8
# Home ranges for full relocation dataset
hr_all <- dat_all %>%
  mutate(n_relocs=map_dbl(data, ~nrow(.))) %>% filter(n_relocs>9) %>%
  mutate(
    hr_mcp = map(data, hr_mcp, levels=hrlevel),
    hr_locoh = map(data, possibly(~ hr_locoh(., n=ceiling(sqrt(nrow(.))), 
                                             levels=hrlevel), NULL)),
    hr_kde = map(data, ~ hr_kde(., trast=make_trast(., factor=1.75), 
                                levels=hrlevel)),
    hr_akde = map(data, ~ hr_akde(., fit_ctmm(., "auto"), levels=hrlevel))
  )
#save(hr_all, file="hr_all.Rdata")
cat(lu(hr_all$id), "fish with all-season home range estimates\n")
## 32 fish with all-season home range estimates
# Home ranges for tracks split by MCP results
hr_season <- dat_season %>%
  mutate(n_relocs=map_dbl(data, ~nrow(.))) %>% filter(n_relocs>9) %>%
  mutate(
    hr_mcp = map(data, hr_mcp, levels=hrlevel),
    hr_locoh = map(data, possibly(~ hr_locoh(., n=ceiling(sqrt(nrow(.))), 
                                             levels=hrlevel), NULL)),
    hr_kde = map(data, ~ hr_kde(., trast=make_trast(., factor=1.75), 
                                levels=hrlevel)),
    hr_akde = map(data, ~ hr_akde(., fit_ctmm(., "auto"), levels=hrlevel))
  )
#save(hr_season, file="hr_season.Rdata")
cat(lu(hr_season$id), "fish with season-specific home range estimates\n")
## 24 fish with season-specific home range estimates
cat(filter(hr_season, season==0) %>% nrow(), "spring home range estimates\n")
## 14 spring home range estimates
cat(filter(hr_season, season==1) %>% nrow(), "summer home range estimates\n")
## 23 summer home range estimates
cat(group_by(hr_season, id) %>% summarize(n=n()) %>% 
      filter(n==2) %>% nrow(), "with both seasons' home range estimates")
## 13 with both seasons' home range estimates

Home range models have now been fit and stored in table objects. Pivot to long format for analysis of home range sizes by estimator.

# pivot to long format
hr_all_1 <- hr_all %>% select(-data) %>% 
  pivot_longer(hr_mcp:hr_akde, names_to="estimator", values_to="hr")
hr_season_1 <- hr_season %>% select(-data) %>% 
  pivot_longer(hr_mcp:hr_akde, names_to="estimator", values_to="hr")
# Functions to calculate areas (in square m)
do_iso <- function(x) hr_isopleths(x) %>% filter(what=="estimate")
clip_iso <- function(x) st_intersection(x, lrwpoly)
# Calculate home range areas
hr_all_1.area <- hr_all_1 %>%
    mutate(hr_area=map(hr, possibly(hr_area, NULL)),
         hr_isopleth=map(hr, possibly(do_iso, NULL)),
         hr_isopleth_clp=map(hr_isopleth, possibly(clip_iso, NULL))
         ) 
hr_season_1.area <- hr_season_1%>% 
  mutate(hr_area=map(hr, possibly(hr_area, NULL)),
         hr_isopleth=map(hr, possibly(do_iso, NULL)),
         hr_isopleth_clp=map(hr_isopleth, possibly(clip_iso, NULL))
         )

Map home ranges

Use the home ranges from aKDE estimation for visualization.

all_akde <- hr_all_1.area %>%  filter(estimator=="hr_akde") %>% 
  select(-hr_isopleth) %>% unnest(hr_isopleth_clp) %>% st_as_sf()
sp_akde <- hr_season_1.area %>%  filter(season==0&estimator=="hr_akde") %>% 
  select(-hr_isopleth) %>% unnest(hr_isopleth_clp) %>% st_as_sf()
su_akde <- hr_season_1.area %>%  filter(season==1&estimator=="hr_akde") %>% 
  select(-hr_isopleth) %>% unnest(hr_isopleth_clp) %>% st_as_sf()
#main plot
plotlocations_f <- ggplot() + 
  geom_tile(data=rw_df %>% filter(!is.na(value)), 
            aes(x=x, y=y), fill='gray90', color='gray90') +
  geom_tile(data=bathy_df %>% filter(!is.na(value)), 
            aes(x=x, y=y, fill=value, color=value)) +
  scale_fill_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
  scale_color_viridis_c(name='Depth (m)', direction=-1, end=0.6) +
    theme_void()
# lines for raw HRs
sp_akde_unclip <- hr_season_1.area %>%  filter(season==0&estimator=="hr_akde") %>% 
  select(-hr_isopleth_clp) %>% unnest(hr_isopleth) %>% st_as_sf() %>% 
  arrange(-area)
su_akde_unclip <- hr_season_1.area %>%  filter(season==1&estimator=="hr_akde") %>% 
  select(-hr_isopleth_clp) %>% unnest(hr_isopleth) %>% st_as_sf() %>% 
  arrange(-area)
# plot fill for clipped but outlines for raw HRs
p_sp_akde_f <- plotlocations_f +
  ggtitle("Spring") +
  geom_sf(data=sp_akde_unclip, alpha=0.2, fill="gold", color=NA) +
  geom_sf(data=sp_akde_unclip, alpha=0.2, fill=NA, color="white") +
  geom_tile(data=lrw_mask_df %>% filter(!is.na(value)), 
            aes(x=x, y=y), fill="white", color="white") +
  geom_sf(data=st_centroid(sp_akde_unclip), shape=3, color="gray20") +
  coord_sf(xlim=st_bbox(olb)[c(1,3)], ylim=st_bbox(olb)[c(2,4)]) 
## Warning: st_centroid assumes attributes are constant over geometries
p_su_akde_f <- plotlocations_f +
  ggtitle("Summer") +
  geom_sf(data=su_akde, alpha=0.2, fill="gold", color=NA) +
  geom_sf(data=su_akde, alpha=0.2, fill=NA, color="white") +
  geom_sf(data=st_centroid(su_akde_unclip), shape=3, color="gray20") +
  coord_sf(xlim=st_bbox(olb)[c(1,3)], ylim=st_bbox(olb)[c(2,4)]) +
    ggspatial::annotation_north_arrow(which_north="true", 
                                      pad_y = unit(0.635, "cm"),
                                      height = unit(1., "cm"),
                                      width = unit(1., "cm"),
        style = ggspatial::north_arrow_fancy_orienteering) +
    ggspatial::annotation_scale(location = "bl", width_hint = 0.5) 
## Warning: st_centroid assumes attributes are constant over geometries
#build plot
combined <- p_sp_akde_f / p_su_akde_f + plot_layout(guides="collect")
ggsave("fig6_pHRS.png", combined, width=6.5, height=4.3, dpi=600, bg="white")
knitr::include_graphics("fig6_pHRS.png")

Plot home range sizes

Plot home range size by estimator and sex for comparison. Home ranges are expressed in meters because that is the basic unit of the coordinate reference system, but here we convert to square km.

Aseasonal home range

# function to derive home range sizes
fci <- function(hr.area){
    hr.area1 <- hr.area %>% 
    unnest(cols=hr_area) %>%
    mutate(area = area / 1e6) 
  hr.area1$est_lab<-factor(hr.area1$estimator, 
                           levels=c("hr_mcp", "hr_locoh", "hr_kde", "hr_akde"), 
                           labels=c("MCP", "LoCoH", "KDE", "aKDE"))
  #remove list column in order to pivot wider and separate estimates from ci
  hr.area2 <- hr.area1 %>% select(-hr) %>%
    pivot_wider(names_from=what, values_from=area)
  #calculate upper and lower confidence intervals on log scale
  ci2 <- hr.area2 %>% group_by(est_lab, sex) %>% 
    summarise(m = mean(log(estimate)), 
              sd = sd(log(estimate)),
              se = sd(log(estimate)) / sqrt(n()), 
              me = qt(0.975, n() - 1) * se, 
              lci = m - me, 
              uci = m + me,
              n=n())
  return(list(hr.area=hr.area2, ci.area=ci2))
}
# function to plot home range sizes
fplot.range <- function(hr.area, ci.area){
  # gj edits to p1
  p1.2 <- hr.area %>% 
    ggplot(aes(est_lab, estimate, shape = sex)) + 
    scale_shape_manual(values=c(2, 5)) +
    geom_pointrange(aes(x = est_lab, y = estimate, ymin = `lci (0.95)`, 
                        ymax = `uci (0.95)`, shape = sex), inherit.aes = FALSE, 
                    position = position_dodge2(width = 0.5)) +
    theme_light() + 
    theme(axis.title.x = element_blank()) +
    labs(y = expression(paste("HRS [", km^2, "]")))
  p2.2 <- hr.area %>% 
    ggplot(aes(est_lab, log(estimate), shape = sex)) + 
    scale_shape_manual(values=c(2, 5)) +
    geom_jitter(alpha = 0.5, position = position_dodge2(width = 0.5)) +
    geom_pointrange(aes(x = est_lab, y = m, ymin = lci, ymax = uci, shape=sex), 
                    data = ci.area, inherit.aes = FALSE, 
                    position = position_dodge2(width = 0.5),
                    size=1) +
    theme_light() + 
    theme(axis.title.x = element_blank()) +
    labs(y = expression(paste("HRS [", km^2, "]")))
  # need confidence interval of the log response ratio log(meanM/meanF)
  # se(log(barm)-log(barf) = sqrt( s^2 /(n*barx)))
  hr.logratio <- hr.area %>% 
    nest(data=-c(est_lab, sex)) %>% 
    mutate(barx=map_dbl(data, ~mean(log(.$estimate))), 
           SD=map_dbl(data, ~sd(log(.$estimate))),
           n=map_dbl(data, ~nrow(.)),
           Sn=(SD^2)/n) %>%
    nest(data=-est_lab) %>%
    mutate(lr=map_dbl(data, ~ .$barx[.$sex=="M"] - .$barx[.$sex=="F"]),
           SElr=map_dbl(data, ~ sqrt(.$Sn[.$sex=="M"] + .$Sn[.$sex=="F"])),
           df=map_dbl(data, ~ .$n[.$sex=="M"]+.$n[.$sex=="F"]-2),
           ucilr=lr+qt(0.975, df)*SElr,
           lcilr=lr-qt(0.975, df)*SElr)
  p3.2 <- hr.logratio %>% 
    ggplot(aes(est_lab, lr)) + 
    geom_pointrange(aes(ymin=lcilr, ymax=ucilr)) + 
    geom_hline(yintercept = 0, lty = 2) +
    theme_light() +
    theme(axis.title.x = element_blank()) +
    labs(y = expression(log(HRS[m]) - log(HRS[f])))
  
  (p1.2 / p2.2 / p3.2) + plot_layout(heights = c(0.5, 0.5, 0.5)) + 
    plot_annotation(tag_levels = "A", tag_suffix = ")")
}
# Run function for non-seasonal home range areas
(hr_all_out<- fci(hr_all_1.area))
## `summarise()` has grouped output by 'est_lab'. You can override using the
## `.groups` argument.
## $hr.area
## # A tibble: 128 × 14
##    id    sex      tl    kg rel       n_relocs estimator level hr_isopleth 
##    <chr> <chr> <dbl> <dbl> <chr>        <dbl> <chr>     <dbl> <list>      
##  1 12    M       545  1.38 Bradburys       24 hr_mcp     0.8  <sf [1 × 4]>
##  2 12    M       545  1.38 Bradburys       24 hr_locoh   0.79 <sf [1 × 4]>
##  3 12    M       545  1.38 Bradburys       24 hr_kde     0.8  <sf [1 × 4]>
##  4 12    M       545  1.38 Bradburys       24 hr_akde    0.8  <sf [1 × 4]>
##  5 14    M       635  1.9  Bradburys       55 hr_mcp     0.8  <sf [1 × 4]>
##  6 14    M       635  1.9  Bradburys       55 hr_locoh   0.85 <sf [1 × 4]>
##  7 14    M       635  1.9  Bradburys       55 hr_kde     0.8  <sf [1 × 4]>
##  8 14    M       635  1.9  Bradburys       55 hr_akde    0.8  <sf [1 × 4]>
##  9 15    F       606  1.8  Bradburys       36 hr_mcp     0.8  <sf [1 × 4]>
## 10 15    F       606  1.8  Bradburys       36 hr_locoh   0.81 <sf [1 × 4]>
## # ℹ 118 more rows
## # ℹ 5 more variables: hr_isopleth_clp <list>, est_lab <fct>, estimate <dbl>,
## #   `lci (0.95)` <dbl>, `uci (0.95)` <dbl>
## 
## $ci.area
## # A tibble: 8 × 9
## # Groups:   est_lab [4]
##   est_lab sex        m    sd    se    me    lci    uci     n
##   <fct>   <chr>  <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl> <int>
## 1 MCP     F      1.27   2.64 0.704 1.52  -0.254  2.79     14
## 2 MCP     M      0.199  1.98 0.466 0.984 -0.785  1.18     18
## 3 LoCoH   F     -0.635  1.71 0.457 0.988 -1.62   0.353    14
## 4 LoCoH   M     -2.53   1.34 0.316 0.667 -3.20  -1.87     18
## 5 KDE     F      4.01   2.04 0.545 1.18   2.83   5.19     14
## 6 KDE     M      3.11   1.64 0.385 0.813  2.30   3.93     18
## 7 aKDE    F      3.74   2.19 0.585 1.26   2.48   5.01     14
## 8 aKDE    M      2.70   1.61 0.380 0.801  1.90   3.50     18

Spring home range

# Run function for seasonal home range areas split by MCP intervals
(hr_season_out_spring <- hr_season_1.area %>% filter(season==0) %>% fci())
## `summarise()` has grouped output by 'est_lab'. You can override using the
## `.groups` argument.
## $hr.area
## # A tibble: 56 × 15
##    id    sex      tl    kg rel       season n_relocs estimator level hr_isopleth
##    <chr> <chr> <dbl> <dbl> <chr>      <dbl>    <dbl> <chr>     <dbl> <list>     
##  1 12    M       545  1.38 Bradburys      0       11 hr_mcp     0.8  <sf>       
##  2 12    M       545  1.38 Bradburys      0       11 hr_locoh   0.82 <sf>       
##  3 12    M       545  1.38 Bradburys      0       11 hr_kde     0.8  <sf>       
##  4 12    M       545  1.38 Bradburys      0       11 hr_akde    0.8  <sf>       
##  5 14    M       635  1.9  Bradburys      0       13 hr_mcp     0.8  <sf>       
##  6 14    M       635  1.9  Bradburys      0       13 hr_locoh   0.77 <sf>       
##  7 14    M       635  1.9  Bradburys      0       13 hr_kde     0.8  <sf>       
##  8 14    M       635  1.9  Bradburys      0       13 hr_akde    0.8  <sf>       
##  9 15    F       606  1.8  Bradburys      0       16 hr_mcp     0.8  <sf>       
## 10 15    F       606  1.8  Bradburys      0       16 hr_locoh   0.81 <sf>       
## # ℹ 46 more rows
## # ℹ 5 more variables: hr_isopleth_clp <list>, est_lab <fct>, estimate <dbl>,
## #   `lci (0.95)` <dbl>, `uci (0.95)` <dbl>
## 
## $ci.area
## # A tibble: 8 × 9
## # Groups:   est_lab [4]
##   est_lab sex        m    sd    se    me     lci     uci     n
##   <fct>   <chr>  <dbl> <dbl> <dbl> <dbl>   <dbl>   <dbl> <int>
## 1 MCP     F     -0.292  2.53 1.03   2.66 -2.95    2.37       6
## 2 MCP     M     -2.19   2.51 0.886  2.09 -4.28   -0.0927     8
## 3 LoCoH   F     -2.22   2.72 1.11   2.86 -5.07    0.641      6
## 4 LoCoH   M     -3.50   2.84 1.00   2.37 -5.87   -1.13       8
## 5 KDE     F      2.44   2.28 0.932  2.40  0.0460  4.84       6
## 6 KDE     M      0.965  1.33 0.469  1.11 -0.144   2.07       8
## 7 aKDE    F      1.84   2.48 1.01   2.60 -0.761   4.44       6
## 8 aKDE    M      0.318  1.50 0.531  1.26 -0.937   1.57       8

Fall home range

# Run function for seasonal home range areas split by MCP intervals
(hr_season_out_su <- hr_season_1.area %>% filter(season==1) %>% fci())
## `summarise()` has grouped output by 'est_lab'. You can override using the
## `.groups` argument.
## $hr.area
## # A tibble: 92 × 15
##    id    sex      tl    kg rel       season n_relocs estimator level hr_isopleth
##    <chr> <chr> <dbl> <dbl> <chr>      <dbl>    <dbl> <chr>     <dbl> <list>     
##  1 29    M       617  2.1  Bradburys      1       26 hr_mcp     0.8  <sf>       
##  2 29    M       617  2.1  Bradburys      1       26 hr_locoh   0.77 <sf>       
##  3 29    M       617  2.1  Bradburys      1       26 hr_kde     0.8  <sf>       
##  4 29    M       617  2.1  Bradburys      1       26 hr_akde    0.8  <sf>       
##  5 12    M       545  1.38 Bradburys      1       13 hr_mcp     0.8  <sf>       
##  6 12    M       545  1.38 Bradburys      1       13 hr_locoh   0.92 <sf>       
##  7 12    M       545  1.38 Bradburys      1       13 hr_kde     0.8  <sf>       
##  8 12    M       545  1.38 Bradburys      1       13 hr_akde    0.8  <sf>       
##  9 15    F       606  1.8  Bradburys      1       11 hr_mcp     0.8  <sf>       
## 10 15    F       606  1.8  Bradburys      1       11 hr_locoh   0.91 <sf>       
## # ℹ 82 more rows
## # ℹ 5 more variables: hr_isopleth_clp <list>, est_lab <fct>, estimate <dbl>,
## #   `lci (0.95)` <dbl>, `uci (0.95)` <dbl>
## 
## $ci.area
## # A tibble: 8 × 9
## # Groups:   est_lab [4]
##   est_lab sex        m    sd    se    me    lci     uci     n
##   <fct>   <chr>  <dbl> <dbl> <dbl> <dbl>  <dbl>   <dbl> <int>
## 1 MCP     F     -0.879  2.81 0.812  1.79 -2.67   0.909     12
## 2 MCP     M     -1.50   2.03 0.613  1.37 -2.87  -0.139     11
## 3 LoCoH   F     -1.73   2.63 0.759  1.67 -3.40  -0.0604    12
## 4 LoCoH   M     -2.86   1.56 0.470  1.05 -3.90  -1.81      11
## 5 KDE     F      1.57   3.12 0.901  1.98 -0.416  3.55      12
## 6 KDE     M      0.590  2.31 0.696  1.55 -0.959  2.14      11
## 7 aKDE    F      1.06   3.00 0.866  1.91 -0.843  2.97      12
## 8 aKDE    M      0.182  2.21 0.667  1.49 -1.31   1.67      11

HRs for appendix

  p1.1 <- hr_season_out_spring$hr.area %>% 
    ggplot(aes(est_lab, estimate, shape = sex)) + 
    scale_shape_manual(values=c(2, 5)) +
    geom_pointrange(aes(x = est_lab, y = estimate, ymin = `lci (0.95)`, 
                        ymax = `uci (0.95)`, shape = sex), inherit.aes = FALSE, 
                    position = position_dodge2(width = 0.5)) +
    theme_light() + 
    theme(axis.title.x = element_blank()) +
    labs(y = expression(paste("HRS [", km^2, "]")))
  p1.2 <- hr_season_out_su$hr.area %>% 
    ggplot(aes(est_lab, estimate, shape = sex)) +
    scale_shape_manual(values=c(2, 5)) +
    geom_pointrange(aes(x = est_lab, y = estimate, ymin = `lci (0.95)`, 
                        ymax = `uci (0.95)`, shape = sex), inherit.aes = FALSE, 
                    position = position_dodge2(width = 0.5)) +
    theme_light() + 
    theme(axis.title.x = element_blank()) +
    labs(y = expression(paste("HRS [", km^2, "]")))
  p2.1 <- hr_season_out_spring$hr.area %>% 
    ggplot(aes(est_lab, log(estimate), shape = sex)) + 
    scale_shape_manual(values=c(2, 5)) +
    geom_jitter(alpha = 0.5, position = position_dodge2(width = 0.5)) +
    geom_pointrange(aes(x = est_lab, y = m, ymin = lci, ymax = uci, shape=sex), 
                    data = hr_season_out_spring$ci.area, inherit.aes = FALSE, 
                    position = position_dodge2(width = 0.5),
                    size=1) +
    theme_light() + 
    theme(axis.title.x = element_blank()) +
    labs(y = expression(paste("HRS [", km^2, "]")))
  p2.2 <- hr_season_out_su$hr.area %>% 
    ggplot(aes(est_lab, log(estimate), shape = sex)) + 
    scale_shape_manual(values=c(2, 5)) +
    geom_jitter(alpha = 0.5, position = position_dodge2(width = 0.5)) +
    geom_pointrange(aes(x = est_lab, y = m, ymin = lci, ymax = uci, shape=sex), 
                    data = hr_season_out_su$ci.area, inherit.aes = FALSE, 
                    position = position_dodge2(width = 0.5),
                    size=1) +
    theme_light() + 
    theme(axis.title.x = element_blank()) +
    labs(y = expression(paste("HRS [", km^2, "]")))
  # need confidence interval of the log response ratio log(meanM/meanF)
  # se(log(barm)-log(barf) = sqrt( s^2 /(n*barx)))
  hr.logratio1 <- hr_season_out_spring$hr.area %>% 
    nest(data=-c(est_lab, sex)) %>% 
    mutate(barx=map_dbl(data, ~mean(log(.$estimate))), 
           SD=map_dbl(data, ~sd(log(.$estimate))),
           n=map_dbl(data, ~nrow(.)),
           Sn=(SD^2)/n) %>%
    nest(data=-est_lab) %>%
    mutate(lr=map_dbl(data, ~ .$barx[.$sex=="M"] - .$barx[.$sex=="F"]),
           SElr=map_dbl(data, ~ sqrt(.$Sn[.$sex=="M"] + .$Sn[.$sex=="F"])),
           df=map_dbl(data, ~.$n[.$sex=="M"]+.$n[.$sex=="F"]-2),
           ucilr=lr+qt(0.975, df)*SElr,
           lcilr=lr-qt(0.975, df)*SElr)
  hr.logratio2 <- hr_season_out_su$hr.area %>% 
    nest(data=-c(est_lab, sex)) %>% 
    mutate(barx=map_dbl(data, ~mean(log(.$estimate))), 
           SD=map_dbl(data, ~sd(log(.$estimate))),
           n=map_dbl(data, ~nrow(.)),
           Sn=(SD^2)/n) %>%
    nest(data=-est_lab) %>%
    mutate(lr=map_dbl(data, ~ .$barx[.$sex=="M"] - .$barx[.$sex=="F"]),
           SElr=map_dbl(data, ~ sqrt(.$Sn[.$sex=="M"] + .$Sn[.$sex=="F"])),
           df=map_dbl(data, ~.$n[.$sex=="M"]+.$n[.$sex=="F"]-2),
           ucilr=lr+qt(0.975, df)*SElr,
           lcilr=lr-qt(0.975, df)*SElr)
  p3.1 <- hr.logratio1 %>% 
    ggplot(aes(est_lab, lr)) + 
    geom_pointrange(aes(ymin=lcilr, ymax=ucilr)) + 
    geom_hline(yintercept = 0, lty = 2) +
    theme_light() +
    theme(axis.title.x = element_blank()) +
    labs(y = expression(log(HRS[m]) - log(HRS[f])))
  p3.2 <- hr.logratio2 %>% 
    ggplot(aes(est_lab, lr)) + 
    geom_pointrange(aes(ymin=lcilr, ymax=ucilr)) + 
    geom_hline(yintercept = 0, lty = 2) +
    theme_light() +
    theme(axis.title.x = element_blank()) +
    labs(y = expression(log(HRS[m]) - log(HRS[f])))
  # build plot
  pspring <- (p1.1 + p2.1 + p3.1) + 
    plot_layout(ncol=1)
  psummer <- (p1.2 + p2.2 + p3.2) + 
    plot_layout(ncol=1)
  pcombo <- (pspring | psummer) + 
    plot_layout(ncol=2) + 
    plot_annotation(tag_levels = "A", tag_suffix = ")")
ggsave("figS4_HRs.png", pcombo, width=9, height=5.3, units="in", dpi=600, 
       bg="white")
## Warning: Removed 42 rows containing missing values (`geom_segment()`).
## Warning: Removed 69 rows containing missing values (`geom_segment()`).
knitr::include_graphics("figS4_HRs.png")

Aggregate and summarize HR estimates

Create a table summarizing home range size esitmates by season, sex, and estimator

(hr.summary.table <- bind_rows(mutate(hr_all_out$hr.area, season=99),
                              hr_season_out_spring$hr.area, 
                              hr_season_out_su$hr.area) %>% 
  mutate(sex=factor(sex), 
         season=factor(season),
         estimator=factor(estimator)) %>% 
  dplyr::select(season, id:level, est_lab, estimate) %>% 
  group_by(sex, estimator, season) %>% 
  summarise(m = mean(log(estimate)), 
            sd = sd(log(estimate)),
            se = sd(log(estimate)) / sqrt(n()), 
            me = qt(0.975, n() - 1) * se, 
            lci = m - me, 
            uci = m + me,
            n=n()) %>% 
  mutate(hat_hr=exp(m), hat_lci=exp(lci), hat_uci=exp(uci)) %>% 
  arrange(estimator, season, sex))
## `summarise()` has grouped output by 'sex', 'estimator'. You can override using
## the `.groups` argument.
## # A tibble: 24 × 13
## # Groups:   sex, estimator [8]
##    sex   estimator season     m    sd    se    me     lci   uci     n hat_hr
##    <fct> <fct>     <fct>  <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl> <int>  <dbl>
##  1 F     hr_akde   0      1.84   2.48 1.01  2.60  -0.761   4.44     6   6.28
##  2 M     hr_akde   0      0.318  1.50 0.531 1.26  -0.937   1.57     8   1.37
##  3 F     hr_akde   1      1.06   3.00 0.866 1.91  -0.843   2.97    12   2.90
##  4 M     hr_akde   1      0.182  2.21 0.667 1.49  -1.31    1.67    11   1.20
##  5 F     hr_akde   99     3.74   2.19 0.585 1.26   2.48    5.01    14  42.1 
##  6 M     hr_akde   99     2.70   1.61 0.380 0.801  1.90    3.50    18  14.9 
##  7 F     hr_kde    0      2.44   2.28 0.932 2.40   0.0460  4.84     6  11.5 
##  8 M     hr_kde    0      0.965  1.33 0.469 1.11  -0.144   2.07     8   2.63
##  9 F     hr_kde    1      1.57   3.12 0.901 1.98  -0.416   3.55    12   4.79
## 10 M     hr_kde    1      0.590  2.31 0.696 1.55  -0.959   2.14    11   1.80
## # ℹ 14 more rows
## # ℹ 2 more variables: hat_lci <dbl>, hat_uci <dbl>
hr.summary.table %>% 
  mutate(cell=paste0(round(hat_hr, 3), "(", 
             round(hat_lci, 3), " - ", 
             round(hat_uci, 3), ")")) %>% 
  dplyr::select(sex, estimator, season, cell) %>% 
  pivot_wider(names_from=c(season, sex), values_from=cell) %>% 
  kableExtra::kable()
estimator 0_F 0_M 1_F 1_M 99_F 99_M
hr_akde 6.282(0.467 - 84.445) 1.375(0.392 - 4.826) 2.897(0.43 - 19.495) 1.199(0.271 - 5.306) 42.121(11.891 - 149.195) 14.938(6.706 - 33.277)
hr_kde 11.499(1.047 - 126.279) 2.626(0.866 - 7.961) 4.787(0.66 - 34.739) 1.805(0.383 - 8.502) 55.238(17.023 - 179.24) 22.509(9.982 - 50.756)
hr_locoh 0.109(0.006 - 1.899) 0.03(0.003 - 0.324) 0.177(0.033 - 0.941) 0.057(0.02 - 0.164) 0.53(0.197 - 1.423) 0.079(0.041 - 0.155)
hr_mcp 0.747(0.052 - 10.672) 0.112(0.014 - 0.911) 0.415(0.069 - 2.481) 0.222(0.057 - 0.871) 3.552(0.776 - 16.266) 1.22(0.456 - 3.265)

Home range size analysis

hr2.area.season <- bind_rows(hr_season_out_spring$hr.area, 
                             hr_season_out_su$hr.area) %>% 
  mutate(sex=factor(sex), 
         season=factor(season),
         estimator=factor(estimator)) %>% 
  dplyr::select(id:level, est_lab, estimate)
m.hr <- lmer(log(estimate) ~ sex*season + (1|est_lab), data=hr2.area.season)
anova(m.hr, ddf="Satterthwaite")
## Type III Analysis of Variance Table with Satterthwaite's method
##            Sum Sq Mean Sq NumDF DenDF F value   Pr(>F)   
## sex        51.361  51.361     1   141  9.0508 0.003111 **
## season      0.471   0.471     1   141  0.0830 0.773637   
## sex:season  3.537   3.537     1   141  0.6233 0.431163   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
margdf <- emmip(m.hr, ~ sex|season, lmer.df="satterthwaite", plotit=FALSE)
pd <- position_dodge(.3)
p <- ggplot(data=margdf, aes(x=season, y=exp(yvar), fill=sex, group=sex)) +
  geom_errorbar(aes(ymin=exp(yvar-SE), ymax=exp(yvar+SE)), 
                color='black', width=.2, position=pd) +
  geom_line(position=pd) +
  geom_point(color='black', size=8, position=pd, shape=21) +
  scale_fill_grey() +
  ylab(expression("Home Range Size "~(km^2))) +
  scale_y_log10() +
  scale_x_discrete(labels=c("0"="Spring", "1"="Summer")) +
  theme_classic() +
  guides(shape=guide_legend(title="Sex")) + 
  theme(axis.title=element_text(size = rel(1.5)),
        axis.text=element_text(size = rel(1.5), color="black"),
        legend.title=element_text(size = rel(1.3)),
        legend.text=element_text(size = rel(1.3)),
        axis.title.x=element_blank())
ggsave("fig7_Ameans.png", p, width=6.5, height=4, units="in", dpi=600, 
       bg="white")
knitr::include_graphics("fig7_Ameans.png")

A tabular summary of regression model output for the appendix.

tab_model(m.hr, df.method="satterthwaite")
  log(estimate)
Predictors Estimates CI p
(Intercept) 0.44 -2.16 – 3.05 0.675
sex [M] -1.54 -2.82 – -0.27 0.018
season [1] -0.44 -1.62 – 0.74 0.463
sex [M] × season [1] 0.64 -0.97 – 2.25 0.431
Random Effects
σ2 5.67
τ00 est_lab 3.00
ICC 0.35
N est_lab 4
Observations 148
Marginal R2 / Conditional R2 0.039 / 0.371

Summer resorts

First, establish a grid over which we may produce count summaries of relocations per fish per season per year. We specified a hexagonal-cell grid (area ~2.5 km per cell) over Oneida Lake and retained all cells in which any fish was observed during our study. Results of seasonal home range analysis suggested that a 2.5 km^2 grid cell could contain the average home range size in spring or summer for males or females.

# create hexagonal grid over Oneida Lake detection locations
g <- hr_season %>% unnest(data) %>% 
  mutate(year=year(t_), yearf=factor(year)) %>% 
  dplyr::select(id, sex, season, year, yearf, x_:t_) %>% 
  st_as_sf(coords=c(y="x_", x="y_"), remove=FALSE, crs=crs_utm18n) %>% 
  st_make_grid(square=FALSE, n = c(10,10)) %>% st_sf() %>% 
  filter(st_intersects(., summarize(fishobs) %>% st_cast("MULTIPOINT"), 
                       sparse=FALSE) %>% as.vector()) %>% 
  mutate(cell=1:nrow(.), area=st_area(.))

Starting from the data table hr_season, create a new table hr_season_summerresort that adds columns for data list cells subset by year (data09, data10, data11) and hr_mcp list cells by year (hr_mcp09, hr_mcp10, hr_mcp11).

hr_summerresort <- hr_season %>% 
  mutate(data09=map(data, ~subset(., year(.$t_)==2009)),
         n09=map_dbl(data09, ~nrow(.)),
         hr_mcp09=ifelse(n09>2, map_dbl(data09, ~nrow(.), levels=1), NA),
         data10=map(data, ~subset(., year(.$t_)==2010)),
         n10=map_dbl(data10, ~nrow(.)),
         hr_mcp10=ifelse(n10>2, map_dbl(data10, ~nrow(.), levels=1), NA),
         data11=map(data, ~subset(., year(.$t_)==2011)),
         n11=map_dbl(data11, ~nrow(.)),
         hr_mcp11=ifelse(n11>2, map_dbl(data11, ~nrow(.), levels=1), NA),
         datacell=map(data, 
                      ~st_intersection(
                        st_as_sf(., coords=c(y="x_", x="y_"), 
                                 remove=FALSE, 
                                 crs=crs_utm18n), 
                        g)
                      ),
         cellct=map(datacell, 
                    ~st_drop_geometry(.) %>%
                      mutate(year=year(t_), 
                             yearf=factor(year, 
                                          levels=c("2009", "2010", "2011"))) %>% 
                      group_by(year, yearf, cell) %>% 
                      summarize(n=n()) %>% 
                      ungroup()
                    ),
         cellp=map(cellct, 
                   ~select(., -yearf) %>% 
                     full_join(expand.grid(year=c(2009:2011),
                                           cell=unique(.$cell)),
                               by=join_by(year, cell)) %>%
                     replace(is.na(.), 0) %>% 
                     group_by(year) %>% 
                     mutate(p=n/sum(n)) %>% ungroup() %>%
                     pivot_wider(id_cols=cell,
                                 names_from=year, 
                                 values_from=p, 
                                 names_prefix='p') %>% 
                     mutate(srpq0910=sqrt(p2009*p2010),
                            srpq1011=sqrt(p2011*p2010)) 
                   ),
         bc0910=map_dbl(cellp, ~summarize(., sum(srpq0910)) %>% pull()),
         bc1011=map_dbl(cellp, ~summarize(., sum(srpq1011)) %>% pull())
         ) %>% 
  rowwise() %>% 
  mutate(bcmean=mean(c(bc0910, bc1011), na.rm=TRUE))
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## `summarise()` has grouped output by 'year', 'yearf'. You can override using the
## `.groups` argument.
## Warning: There were 37 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `datacell = map(...)`.
## Caused by warning:
## ! attribute variables are assumed to be spatially constant throughout all geometries
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 36 remaining warnings.

Plot summer locations and cell distributions for individual fish. Ladies first.

scale_color_years <- function(...){
    ggplot2:::manual_scale(
        'color', 
        values = setNames(c("#f98e09", "#57106e", "#000004"), 
                          #c("2009", "2010", "2011")),
                          c("2011", "2010", "2009")), 
        ...
    )
}
hr_summerresortM <- hr_summerresort %>% 
  filter(sex=="M", season==1&!is.na(bcmean)) %>% arrange(-bcmean)
psrlistM <- lapply(1:nrow(hr_summerresortM), function(i){
  sra <- hr_summerresortM[i,] %>% 
    unnest(data) %>% 
    st_as_sf(coords=c("x_", "y_"), crs=crs_utm18n)
  this <-  ggplot(sra) +
    geom_text(aes(x=411000, y=4790200, label= 
              paste0("ID = ", hr_summerresortM[i,"id"],
                     ", BC = ", round(hr_summerresortM[i,"bcmean"], 2))),
              size=2.5, hjust=0) +
    geom_sf(data=g, fill=NA, color="grey90") +
    geom_text(aes(x=408000, y=4782000), 
             label=paste0("2009, n=", sum(year(sra$t_)==2009)),
             size=2, hjust=0,  parse=FALSE, color="#000004") +
    geom_text(aes(x=408000, y=4781000), 
             label=paste0("2010, n=", sum(year(sra$t_)==2010)),
             size=2, hjust=0,  parse=FALSE, color="#57106e") +
    geom_text(aes(x=408000, y=4780000), 
             label=paste0("2011, n=", sum(year(sra$t_)==2011)),
             size=2, hjust=0,  parse=FALSE, color="#f98e09") +
    geom_sf(data=olb, fill=NA, color="grey50") +
    geom_sf(aes(color=factor(year(t_))), shape=3, size=3, alpha=1) + 
    scale_color_years(name="Year") +
    coord_sf(xlim=st_bbox(g)[c(1,3)], ylim=st_bbox(g)[c(2,4)]) +
    theme_void() + 
    theme(plot.margin=unit(c(0,0,0,0), "cm"),
          plot.title=element_text(size=10)) + 
    theme(legend.position="none")
  return(this)
})
for (x in 1:(length(psrlistM)-1)){
psrlistM[[x]] <- psrlistM[[x]] + 
    guides(shape="none", color="none", fill="none") + 
    theme(legend.position="none")
}
fig9 <- wrap_plots(psrlistM) + plot_layout(ncol=3, guides="collect") + 
  plot_annotation(tag_levels="a", tag_suffix=")") 
ggsave("fig9_pBCM.png", fig9, width=6.5, height=6, dpi=300, bg="white")
knitr::include_graphics("fig9_pBCM.png")

hr_summerresortF <- hr_summerresort %>% 
  filter(sex=="F", season==1&!is.na(bcmean)) %>% arrange(-bcmean)
psrlistF <- lapply(1:nrow(hr_summerresortF), function(i){
  sra <- hr_summerresortF[i,] %>% 
    unnest(data) %>% 
    st_as_sf(coords=c("x_", "y_"), crs=crs_utm18n)
  this <-  ggplot(sra) +
    geom_text(aes(x=411000, y=4790200, label= 
              paste0("ID = ", hr_summerresortF[i,"id"],
                     ", BC = ", round(hr_summerresortF[i,"bcmean"], 2))),
              size=2.5, hjust=0) +
    geom_sf(data=g, fill=NA, color="grey90") +
    geom_text(aes(x=408000, y=4782000), 
             label=paste0("2009, n=", sum(year(sra$t_)==2009)),
             size=2, hjust=0,  parse=FALSE, color="#000004") +
    geom_text(aes(x=408000, y=4781000), 
             label=paste0("2010, n=", sum(year(sra$t_)==2010)),
             size=2, hjust=0,  parse=FALSE, color="#57106e") +
    geom_text(aes(x=408000, y=4780000), 
             label=paste0("2011, n=", sum(year(sra$t_)==2011)),
             size=2, hjust=0,  parse=FALSE, color="#f98e09") +
    geom_sf(data=olb, fill=NA, color="grey50") +
    geom_sf(aes(color=factor(year(t_))), shape=3, size=3, alpha=1) + 
    scale_color_years(name="Year") +
    coord_sf(xlim=st_bbox(g)[c(1,3)], ylim=st_bbox(g)[c(2,4)]) +
    theme_void() + 
    theme(plot.margin=unit(c(0,0,0,0), "cm"),
          plot.title=element_text(size=10)) + 
    theme(legend.position="none")
  return(this)
})
for (x in 1:(length(psrlistF)-1)){
psrlistF[[x]] <- psrlistF[[x]] + 
    guides(shape="none", color="none", fill="none") + 
    theme(legend.position="none")
}
fig10 <- wrap_plots(psrlistF) + plot_layout(ncol=3, guides="collect") + 
  plot_annotation(tag_levels="a", tag_suffix=")") 
ggsave("fig10_BCF.png", fig10, width=6.5, height=6, dpi=300, bg="white")
knitr::include_graphics("fig10_BCF.png")

fig8 <- hr_summerresort %>% filter(season==1&!is.na(bcmean)) %>% 
  ggplot(aes(x=sex, y=bcmean)) + 
  geom_point(aes(fill=sex), size=3, shape=21, color="grey20",
             position=position_jitter(width=0.2)) + 
  geom_boxplot(outlier.color=NA, fill=NA, color="grey20") +
  scale_fill_grey(name="Sex") +
  ylab("Bhattacharyya's coeffcient (BC)") + xlab("") +
  ylim(0,1) +
  theme_classic()
ggsave("fig8_pba_overlap.png", fig8, width=3, height=3, dpi=300, bg="white")
## Warning: Removed 1 rows containing missing values (`geom_point()`).
knitr::include_graphics("fig8_pba_overlap.png")

hr_summerresort %>% filter(season==1&!is.na(bcmean)) %>%
  group_by(sex) %>% 
  summarize(n=n(), mean=mean(bcmean), sd=sd(bcmean), 
            q25=quantile(bcmean, probs=c(0.25)), 
            q75=quantile(bcmean, probs=c(0.75)), 
            median=median(bcmean), min=min(bcmean), max=max(bcmean), )
## # A tibble: 2 × 9
##   sex       n  mean    sd   q25   q75 median   min   max
##   <chr> <int> <dbl> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>
## 1 F        12 0.672 0.294 0.498 0.895  0.752 0.167 1    
## 2 M        11 0.693 0.360 0.666 0.939  0.822 0     0.964
Sys.time()-t0
## Time difference of 6.737476 mins
testjags()
## You are using R version 4.3.1 (2023-06-16 ucrt) on a windows machine,
## with the RTerm GUI
## JAGS version 4.3.1 found successfully using the command '/Program
## Files/JAGS/JAGS-4.3.1/x64/bin/jags-terminal.exe'
## The rjags package is installed
sessionInfo()
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19045)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] DHARMa_0.4.6          sjlabelled_1.2.0      sjmisc_2.8.9         
##  [4] sjPlot_2.8.15         emmeans_1.9.0         lmerTest_3.1-3       
##  [7] lme4_1.1-35.1         Matrix_1.6-0          patchwork_1.2.0      
## [10] runjags_2.2.2-1.1     amt_0.2.1.0           raster_3.6-26        
## [13] sp_2.1-2              sf_1.0-14             readxl_1.4.3         
## [16] kableExtra_1.3.4.9000 ggpubr_0.6.0          ggExtra_0.10.1       
## [19] lubridate_1.9.3       forcats_1.0.0         stringr_1.5.1        
## [22] dplyr_1.1.2           purrr_1.0.2           readr_2.1.4          
## [25] tidyr_1.3.0           tibble_3.2.1          ggplot2_3.4.4        
## [28] tidyverse_2.0.0      
## 
## loaded via a namespace (and not attached):
##   [1] rstudioapi_0.15.0   jsonlite_1.8.8      datawizard_0.9.1   
##   [4] magrittr_2.0.3      ggbeeswarm_0.7.2    TH.data_1.1-2      
##   [7] estimability_1.4.1  farver_2.1.1        nloptr_2.0.3       
##  [10] rmarkdown_2.25      ragg_1.2.7          vctrs_0.6.3        
##  [13] minqa_1.2.6         effectsize_0.8.6    terra_1.7-65       
##  [16] rstatix_0.7.2       webshot_0.5.5       htmltools_0.5.7    
##  [19] broom_1.0.5         cellranger_1.1.0    sass_0.4.8         
##  [22] KernSmooth_2.23-22  bslib_0.6.1         sandwich_3.1-0     
##  [25] zoo_1.8-12          cachem_1.0.8        mime_0.12          
##  [28] lifecycle_1.0.4     pkgconfig_2.0.3     R6_2.5.1           
##  [31] fastmap_1.1.1       rbibutils_2.2.16    shiny_1.8.0        
##  [34] digest_0.6.34       numDeriv_2016.8-1.1 colorspace_2.1-0   
##  [37] RSpectra_0.16-1     textshaping_0.3.7   labeling_0.4.3     
##  [40] fansi_1.0.4         timechange_0.2.0    httr_1.4.7         
##  [43] abind_1.4-5         compiler_4.3.1      proxy_0.4-27       
##  [46] bit64_4.0.5         withr_3.0.0         backports_1.4.1    
##  [49] carData_3.0-5       DBI_1.2.1           performance_0.10.8 
##  [52] highr_0.10          ggsignif_0.6.4      MASS_7.3-60.0.1    
##  [55] sjstats_0.18.2      classInt_0.4-10     loo_2.6.0          
##  [58] tools_4.3.1         units_0.8-4         vipor_0.4.7        
##  [61] lmtest_0.9-40       beeswarm_0.4.0      httpuv_1.6.13      
##  [64] glue_1.6.2          nlme_3.1-162        promises_1.2.1     
##  [67] grid_4.3.1          checkmate_2.3.1     generics_0.1.3     
##  [70] gtable_0.3.4        tzdb_0.4.0          class_7.3-22       
##  [73] hms_1.1.3           xml2_1.3.6          car_3.1-2          
##  [76] utf8_1.2.3          pillar_1.9.0        ctmm_1.2.0         
##  [79] vroom_1.6.3         ggspatial_1.1.9     later_1.3.2        
##  [82] robustbase_0.99-1   splines_4.3.1       lattice_0.22-5     
##  [85] FNN_1.1.4           bit_4.0.5           survival_3.5-7     
##  [88] tidyselect_1.2.0    Gmedian_1.2.7       sfheaders_0.4.3    
##  [91] miniUI_0.1.1.1      knitr_1.45          svglite_2.1.3      
##  [94] xfun_0.39           expm_0.999-9        matrixStats_1.2.0  
##  [97] DEoptimR_1.1-3      stringi_1.7.12      yaml_2.3.8         
## [100] boot_1.3-28.1       evaluate_0.23       codetools_0.2-19   
## [103] cli_3.6.1           parameters_0.21.3   xtable_1.8-4       
## [106] systemfonts_1.0.5   Rdpack_2.6          munsell_0.5.0      
## [109] jquerylib_0.1.4     modelr_0.1.11       Rcpp_1.0.11        
## [112] ggeffects_1.3.4     png_0.1-8           coda_0.19-4        
## [115] parallel_4.3.1      ellipsis_0.3.2      mcp_0.3.3          
## [118] bayestestR_0.13.1   viridisLite_0.4.2   mvtnorm_1.2-4      
## [121] scales_1.3.0        e1071_1.7-13        crayon_1.5.2       
## [124] insight_0.19.7      rlang_1.1.1         rvest_1.0.3        
## [127] multcomp_1.4-25